Kodiak Chum Baseline - 2018 Update

knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)
library(lubridate)
library(leaflet)
library(DT)
library(genepop)
library(abind)
bbind <- function(...) { abind(..., along = 3) }
Sys.setenv(RGL_USE_NULL = TRUE)
library(rgl)
.username <- "krshedd"
.password <- ""
source("~/../R/Functions.GCL.R")

Introduction

Purpose

This is a quick baseline update to Andy Barclay’s 2015 Kodiak chum baseline. This baseline update includes 3 new collections sampled by ADF&G staff (Birch Foster) in 2017:

1. CMBARL17 - Barling Bay Creek - 9/1/2017
2. CMKIAV17 - Kiavak Portage - 9/1/2017
3. CMNATA17 - Natalia Bay Creek - 9/1/2017

Background

Andy’s work in late 2015 was done to assess population genetic structure of chum salmon in SE Kodiak, specifically in the Three Saints Bay area. Kodiak Regional Aquaculture Association (KRAA) wanted to establish a new hatchery in Old Harbor to make early-run chum salmon for a remote release site in Three Saints Bay. KRAA intended to use the Sturgeon broodstock, that they already use for their Kitoi Bay facility. ADF&G regional staff indicated that there are known early-run chum stocks in the Three Saints Bay area, which would present a genetic risk to using Sturgeon River broodstock in the area (similar run timing, but potentially genetically different). Regional staff indicated that Barling Bay Creek in particular had early-run chum salmon. Thus, regional staff collected temporal samples throughout the season from Barling Bay Creek and other nearby systems to get a better understanding of population genetic structure in SE Kodiak to determine whether early-run Barling Bay Creek chum were sufficiently distinct from Sturgeon River chum to be an unaccpetable genetic risk. Ultimately, Andy’s work did show that Barling Bay Creek has significant temporal genetic structure throughout the season (July vs. August vs. September). These results resulted in a memo with a department recommendation to use local broodstock (early-run Barling Bay Creek) rather than Sturgeon River broodstock from Kitoi Bay for the Three Saints Bay project.

This 2018 baseline update will incorporate three additional collections made by ADF&G staff. The results from this baseline update are not going to change the decision regarding use of local broodsource, they are intended to provide some additional pre-hatchery baseline.

Outline

This R Notebook will perform the following tasks:

• Read in genotype data
  • Create LocusControl
  • Read in all project genotypes
• Map of all baseline collections
• Pooling Barling 2015
  • Split CMBARLR15 into early, middle, and late collections (as in the 2015 baseline)
• Determine FailureRate
• Data QA
  • Removing fish missing >= 20% of genotypes
  • Remove within collection duplicates
  • Save final, post-QA genotypes
• Check HWE - collections
• Pooling temporal collections
• Check HWE - populations
• Map of all baseline populations
• Explore genetic structure
  • Generate allele frequency plots
  • Fst tree
  • MDS using Fst

Read in genotype data

Need to create LocusControl for the markers we are using and then read in all the genotype data for the same collections that Andy analyzed in 2015 + the three new collections from 2017 (the update).

Markerset

This analysis is using the 96 SNPs from WASSIP.

CreateLocusControl.GCL(markersuite = "ChumGolden2011_96SNPs", username = .username, password = .password)

[1] "LocusControl created for markersuite='ChumGolden2011_96SNPs'"

loci96 <- LocusControl$locusnames
dir.create("../2018/Objects")
save_objects(c("LocusControl", "loci96"), "../2018/Objects")

Collections

We want to include the same 86 collections used by Andy in his 2015 baseline + the 3 new collections from 2017.

(collections_89 <- read_csv("../2018/collections89.csv") %>% 
  pull("silly"))

Parsed with column specification:
cols(
  silly = col_character()
)

 [1] "CHBAL92UW"    "CHLIJ92UW"    "CMLIJ09UW"    "CMBEL92UW"    "CMCAN09UW"    "CMCAN92UW"    "CMCHI09UW"    "CMCHI96UW"   
 [9] "CMCOL96UW"    "CMDEL96UW"    "CMFOST09UW"   "CMRUB96UW"    "CMRUS09UW"    "CMRUS92UW"    "CMRUS93UW"    "CMSANC09UW"  
[17] "CMSANC96UW"   "CMSMOKH10"    "CMSTE09UW"    "CMSTE92UW"    "CMSTEPR09UW"  "CMWESM93UW"   "CMVOL09UW"    "CMVOL92UW"   
[25] "CMVOL96UW"    "CMZAC92UW"    "CMWESN93UW"   "CMAMBM09UW"   "CMCHIGK09UW"  "CMWESI93UW"   "CMIVAN09UW"   "CMWESL93UW"  
[33] "CMWESK93UW"   "CMKUJUNF09UW" "CMNAKIL08UW"  "CMNECR08UW"   "CMOCEB09UW"   "CMPORTC08UW"  "CMWESJ93UW"   "CMRUDY10"    
[41] "CMWESE93UW"   "CMWESD93UW"   "CMBEARBC09UW" "CMWESH93UW"   "CMBIGRI09UW"  "CMDRYBR09UW"  "CMWESF93UW"   "CMKIAL09UW"  
[49] "CMPASS09UW"   "CMWESG93UW"   "CHAME92UW"    "CHDOG92UW"    "CHKIZ92UW"    "CMKIZH09UW"   "CHUGA92UW"    "CMUGAN09UW"  
[57] "CMBARL09UW"   "CMBARLR15"    "CMBSU92UW"    "CMBSUK09UW"   "CMCOXC09UW"   "CMDEAD09UW"   "CMDOGSC10"    "CMEAGLH09UW" 
[65] "CMWESP93UW"   "CMGULLC09UW"  "CMKARL09UW"   "CMWESA93UW"   "CMKITB09UW"   "CMPAUL15"     "CMMIDWAY15"   "CMHIDDB15"   
[73] "CMKAIUG15"    "CMOLDS10"     "CMPORTNE09UW" "CMROUGH09UW"  "CMRUSSI07UW"  "CMRUSSI09UW"  "CMSITKI09UW"  "CMSPIRU09UW" 
[81] "CMSTU08UW"    "CMSTU09UW"    "CMSTU92"      "CMWKIL08UW"   "CMWKIL09UW"   "CMZACR09UW"   "CMBARL17"     "CMKIAV17"    
[89] "CMNATA17"    

save_objects("collections_89", "../2018/Objects/")

Read in genotype data from LOKI.

LOKI2R.GCL(sillyvec = collections_89, username = .username, password = .password)

CHBAL92UW.gcl created 1 of 89 completed.
CHLIJ92UW.gcl created 2 of 89 completed.
CMLIJ09UW.gcl created 3 of 89 completed.
CMBEL92UW.gcl created 4 of 89 completed.
CMCAN09UW.gcl created 5 of 89 completed.
CMCAN92UW.gcl created 6 of 89 completed.
CMCHI09UW.gcl created 7 of 89 completed.
CMCHI96UW.gcl created 8 of 89 completed.
CMCOL96UW.gcl created 9 of 89 completed.
CMDEL96UW.gcl created 10 of 89 completed.
CMFOST09UW.gcl created 11 of 89 completed.
CMRUB96UW.gcl created 12 of 89 completed.
CMRUS09UW.gcl created 13 of 89 completed.
CMRUS92UW.gcl created 14 of 89 completed.
CMRUS93UW.gcl created 15 of 89 completed.
CMSANC09UW.gcl created 16 of 89 completed.
CMSANC96UW.gcl created 17 of 89 completed.
CMSMOKH10.gcl created 18 of 89 completed.
CMSTE09UW.gcl created 19 of 89 completed.
CMSTE92UW.gcl created 20 of 89 completed.
CMSTEPR09UW.gcl created 21 of 89 completed.
CMWESM93UW.gcl created 22 of 89 completed.
CMVOL09UW.gcl created 23 of 89 completed.
CMVOL92UW.gcl created 24 of 89 completed.
CMVOL96UW.gcl created 25 of 89 completed.
CMZAC92UW.gcl created 26 of 89 completed.
CMWESN93UW.gcl created 27 of 89 completed.
CMAMBM09UW.gcl created 28 of 89 completed.
CMCHIGK09UW.gcl created 29 of 89 completed.
CMWESI93UW.gcl created 30 of 89 completed.
CMIVAN09UW.gcl created 31 of 89 completed.
CMWESL93UW.gcl created 32 of 89 completed.
CMWESK93UW.gcl created 33 of 89 completed.
CMKUJUNF09UW.gcl created 34 of 89 completed.
CMNAKIL08UW.gcl created 35 of 89 completed.
CMNECR08UW.gcl created 36 of 89 completed.
CMOCEB09UW.gcl created 37 of 89 completed.
CMPORTC08UW.gcl created 38 of 89 completed.
CMWESJ93UW.gcl created 39 of 89 completed.
CMRUDY10.gcl created 40 of 89 completed.
CMWESE93UW.gcl created 41 of 89 completed.
CMWESD93UW.gcl created 42 of 89 completed.
CMBEARBC09UW.gcl created 43 of 89 completed.
CMWESH93UW.gcl created 44 of 89 completed.
CMBIGRI09UW.gcl created 45 of 89 completed.
CMDRYBR09UW.gcl created 46 of 89 completed.
CMWESF93UW.gcl created 47 of 89 completed.
CMKIAL09UW.gcl created 48 of 89 completed.
CMPASS09UW.gcl created 49 of 89 completed.
CMWESG93UW.gcl created 50 of 89 completed.
CHAME92UW.gcl created 51 of 89 completed.
CHDOG92UW.gcl created 52 of 89 completed.
CHKIZ92UW.gcl created 53 of 89 completed.
CMKIZH09UW.gcl created 54 of 89 completed.
CHUGA92UW.gcl created 55 of 89 completed.
CMUGAN09UW.gcl created 56 of 89 completed.
CMBARL09UW.gcl created 57 of 89 completed.
CMBARLR15.gcl created 58 of 89 completed.
CMBSU92UW.gcl created 59 of 89 completed.
CMBSUK09UW.gcl created 60 of 89 completed.
CMCOXC09UW.gcl created 61 of 89 completed.
CMDEAD09UW.gcl created 62 of 89 completed.
CMDOGSC10.gcl created 63 of 89 completed.
CMEAGLH09UW.gcl created 64 of 89 completed.
CMWESP93UW.gcl created 65 of 89 completed.
CMGULLC09UW.gcl created 66 of 89 completed.
CMKARL09UW.gcl created 67 of 89 completed.
CMWESA93UW.gcl created 68 of 89 completed.
CMKITB09UW.gcl created 69 of 89 completed.
CMPAUL15.gcl created 70 of 89 completed.
CMMIDWAY15.gcl created 71 of 89 completed.
CMHIDDB15.gcl created 72 of 89 completed.
CMKAIUG15.gcl created 73 of 89 completed.
CMOLDS10.gcl created 74 of 89 completed.
CMPORTNE09UW.gcl created 75 of 89 completed.
CMROUGH09UW.gcl created 76 of 89 completed.
CMRUSSI07UW.gcl created 77 of 89 completed.
CMRUSSI09UW.gcl created 78 of 89 completed.
CMSITKI09UW.gcl created 79 of 89 completed.
CMSPIRU09UW.gcl created 80 of 89 completed.
CMSTU08UW.gcl created 81 of 89 completed.
CMSTU09UW.gcl created 82 of 89 completed.
CMSTU92.gcl created 83 of 89 completed.
CMWKIL08UW.gcl created 84 of 89 completed.
CMWKIL09UW.gcl created 85 of 89 completed.
CMZACR09UW.gcl created 86 of 89 completed.
CMBARL17.gcl created 87 of 89 completed.
CMKIAV17.gcl created 88 of 89 completed.
CMNATA17.gcl created 89 of 89 completed.

Time difference of 2.289525 mins
Time difference of 2.289525 mins

rm(.username, .password)
dir.create("../2018/Genotypes")
dir.create("../2018/Genotypes/original")
save_sillys(collections_89, "../2018/Genotypes/original")

Map of all baseline collections

Pull the latitude/longitude data for collections_89 from OceanAK, need to remove “UW” from all sillys for filtering to get the report.

writeClipboard(paste(str_replace(collections_89, "UW", ""), collapse = ";"))
dir.create("../2018/OceanAK")

Create a map of the baseline collections

read_csv("../2018/OceanAK/collections_89_Just the Lat_Longs.csv") %>% 
  leaflet(width = "100%") %>% 
  addTiles() %>% 
  addMarkers(~ Longitude, ~ Latitude, label = ~ `Silly Code`, clusterOptions = markerClusterOptions())

Parsed with column specification:
cols(
  `Collection ID` = col_integer(),
  `Silly Code` = col_character(),
  Latitude = col_double(),
  Longitude = col_double(),
  `Common Name` = col_character()
)

Create an awesome map of the baseline collections

collections_89_info <- read_csv("../2018/OceanAK/collections_89_Collection Info_mod.csv") 

Parsed with column specification:
cols(
  `Silly Code` = col_character(),
  Region = col_character(),
  Quadrant = col_character(),
  Location = col_character(),
  Latitude = col_double(),
  Longitude = col_double(),
  `Collection Date` = col_character(),
  `Date Received` = col_character()
)

collections_89_info %>% 
  count(Quadrant)

collections_89_info <- collections_89_info %>% 
  mutate(color = case_when(Quadrant == "Peninsula - South" ~ "red",
                           Quadrant == "Chignik" ~ "orange",
                           Quadrant == "Kodiak Mainland" ~ "blue",
                           Quadrant == "Kodiak/Afognak Islands" ~ "purple"))
icons <- awesomeIcons(icon = 'ios-close', iconColor = 'black', library = 'ion', markerColor = collections_89_info$color)
collections_89_info %>% 
  leaflet(width = "100%") %>% 
  addTiles() %>% 
  addAwesomeMarkers(~ Longitude, ~ Latitude, icon = icons, label = ~ `Silly Code`)

Pooling Barling 2015

We need to split CMBARLR15 into early, middle, and late collections (as in the 2015 baseline), as this silly represents 3 distinct collections.

table(CMBARLR15.gcl$attributes$CAPTURE_DATE)

2015-07-09 2015-08-03 2015-09-11 
       104        104         88 

Use AttributesToIDs.GCL to split by date and pool into new silly’s.

# early
CMBARLR15_early_ids <- AttributesToIDs.GCL(silly = "CMBARLR15", attribute = "CAPTURE_DATE", matching = "2015-07-09")
PoolCollections.GCL(collections = "CMBARLR15", loci = loci96, IDs = list(CMBARLR15 = CMBARLR15_early_ids), newname = "CMBARLR15E")

NULL

# middle
CMBARLR15_middle_ids <- AttributesToIDs.GCL(silly = "CMBARLR15", attribute = "CAPTURE_DATE", matching = "2015-08-03")
PoolCollections.GCL(collections = "CMBARLR15", loci = loci96, IDs = list(CMBARLR15 = CMBARLR15_middle_ids), newname = "CMBARLR15M")

NULL

# late
CMBARLR15_late_ids <- AttributesToIDs.GCL(silly = "CMBARLR15", attribute = "CAPTURE_DATE", matching = "2015-09-11")
PoolCollections.GCL(collections = "CMBARLR15", loci = loci96, IDs = list(CMBARLR15 = CMBARLR15_late_ids), newname = "CMBARLR15L")

NULL

save_sillys(c("CMBARLR15E", "CMBARLR15M", "CMBARLR15L"), "../2018/Genotypes/original/")

Now create a new sillyvec with the three temporal 2015 Barling collections split out and the original CMBARLR15 removed. Planning to sort for ease of use.

(collections_91 <- sort(c(collections_89[collections_89 != "CMBARLR15"], "CMBARLR15E", "CMBARLR15M", "CMBARLR15L")))

 [1] "CHAME92UW"    "CHBAL92UW"    "CHDOG92UW"    "CHKIZ92UW"    "CHLIJ92UW"    "CHUGA92UW"    "CMAMBM09UW"   "CMBARL09UW"  
 [9] "CMBARL17"     "CMBARLR15E"   "CMBARLR15L"   "CMBARLR15M"   "CMBEARBC09UW" "CMBEL92UW"    "CMBIGRI09UW"  "CMBSU92UW"   
[17] "CMBSUK09UW"   "CMCAN09UW"    "CMCAN92UW"    "CMCHI09UW"    "CMCHI96UW"    "CMCHIGK09UW"  "CMCOL96UW"    "CMCOXC09UW"  
[25] "CMDEAD09UW"   "CMDEL96UW"    "CMDOGSC10"    "CMDRYBR09UW"  "CMEAGLH09UW"  "CMFOST09UW"   "CMGULLC09UW"  "CMHIDDB15"   
[33] "CMIVAN09UW"   "CMKAIUG15"    "CMKARL09UW"   "CMKIAL09UW"   "CMKIAV17"     "CMKITB09UW"   "CMKIZH09UW"   "CMKUJUNF09UW"
[41] "CMLIJ09UW"    "CMMIDWAY15"   "CMNAKIL08UW"  "CMNATA17"     "CMNECR08UW"   "CMOCEB09UW"   "CMOLDS10"     "CMPASS09UW"  
[49] "CMPAUL15"     "CMPORTC08UW"  "CMPORTNE09UW" "CMROUGH09UW"  "CMRUB96UW"    "CMRUDY10"     "CMRUS09UW"    "CMRUS92UW"   
[57] "CMRUS93UW"    "CMRUSSI07UW"  "CMRUSSI09UW"  "CMSANC09UW"   "CMSANC96UW"   "CMSITKI09UW"  "CMSMOKH10"    "CMSPIRU09UW" 
[65] "CMSTE09UW"    "CMSTE92UW"    "CMSTEPR09UW"  "CMSTU08UW"    "CMSTU09UW"    "CMSTU92"      "CMUGAN09UW"   "CMVOL09UW"   
[73] "CMVOL92UW"    "CMVOL96UW"    "CMWESA93UW"   "CMWESD93UW"   "CMWESE93UW"   "CMWESF93UW"   "CMWESG93UW"   "CMWESH93UW"  
[81] "CMWESI93UW"   "CMWESJ93UW"   "CMWESK93UW"   "CMWESL93UW"   "CMWESM93UW"   "CMWESN93UW"   "CMWESP93UW"   "CMWKIL08UW"  
[89] "CMWKIL09UW"   "CMZAC92UW"    "CMZACR09UW"  

save_objects("collections_91", "../2018/Objects/")

Save sillys post split.

dir.create("../2018/Genotypes/original_split")
save_sillys(collections_91, "../2018/Genotypes/original_split")

Determine FailureRate

project <- "westward_chum_baseline"
loci <- loci96
(failure_rate <- FailureRate.GCL(sillyvec = collections_91))

$`silly_failure_rate`

$locus_failure_rate

$plate_failure_rate

$overall_failure_rate

$plot_silly_failure_rate

$plot_plate_failure_rate

failure_rate_noplots <- failure_rate[1:4]
save_objects("failure_rate_noplots", "../2018/Objects/")

Also calculate sample size by locus and save.

(sample_size_by_locus <- SampSizeByLocus.GCL(sillyvec = collections_91, loci = loci96) %>% 
  as_tibble(rownames = "silly"))

dir.create("../2018/Tables")
write_csv(sample_size_by_locus, "../2018/Tables/sample_size_by_locus.csv")

Data QA

Perform standard data QA processes to filter out untrustworth genotypes.

sample_size_qa <- tibble(silly = collections_91) %>% 
  mutate(genotyped = sapply(silly, function(x) get(paste0(x, ".gcl"))$n))

Missing

Remove individuals missing >=20% of genotypes (i.e. the 80% rule).

miss_loci <- RemoveIndMissLoci.GCL(sillyvec = collections_91, proportion = 0.8)
save_objects("miss_loci", "../2018/Objects/")
# show individuals removed
miss_loci[miss_loci != "None"]

$`CHBAL92UW`
[1] "10" "22"

$CHKIZ92UW
[1] "48" "54" "64" "81"

$CHUGA92UW
[1] "60"

$CMBARL09UW
[1] "92" "95"

$CMBARL17
[1] "63"

$CMBARLR15E
[1] "50" "56" "59" "61" "62" "77"

$CMBARLR15L
 [1] "3"  "16" "17" "19" "20" "42" "46" "61" "76" "85"

$CMBARLR15M
[1] "26" "28"

$CMBEARBC09UW
[1] "66" "84"

$CMCAN09UW
[1] "58" "64" "71"

$CMCAN92UW
[1] "93"

$CMCHI09UW
[1] "16" "52" "95"

$CMCHI96UW
[1] "19" "29" "54"

$CMCHIGK09UW
[1] "38"

$CMCOXC09UW
[1] "39"

$CMFOST09UW
[1] "25" "43" "72" "80" "86" "94"

$CMGULLC09UW
[1] "86"

$CMKAIUG15
[1] "66" "80"

$CMKARL09UW
[1] "15" "25" "45" "58" "60" "66" "91"

$CMKITB09UW
  [1] "96"  "97"  "98"  "99"  "100" "101" "102" "103" "104" "105" "106" "107" "108" "109" "110" "111" "112" "113" "114" "115"
 [21] "116" "117" "118" "119" "120" "121" "122" "123" "124" "125" "126" "127" "128" "129" "130" "131" "132" "133" "134" "135"
 [41] "136" "137" "138" "139" "140" "141" "142" "143" "144" "145" "146" "147" "148" "149" "150" "151" "152" "153" "154" "155"
 [61] "156" "157" "158" "159" "160" "161" "162" "163" "164" "165" "166" "167" "168" "169" "170" "171" "172" "173" "174" "175"
 [81] "176" "177" "178" "179" "180" "181" "182" "183" "184" "185" "186" "187" "188" "189" "190" "191" "192" "193" "194" "195"
[101] "196" "197" "198" "199" "200" "201" "202"

$CMKIZH09UW
[1] "10" "14" "31" "55"

$CMKUJUNF09UW
[1] "14" "15"

$CMLIJ09UW
[1] "8"  "33" "34"

$CMMIDWAY15
[1] "1"

$CMNAKIL08UW
 [1] "96"  "97"  "98"  "99"  "100" "101" "102" "103" "104" "105" "106" "107"

$CMNECR08UW
 [1] "13"  "96"  "97"  "98"  "99"  "100" "101" "102" "103" "104" "105" "106" "107" "108" "109" "110" "111" "112"

$CMPASS09UW
[1] "41"

$CMPAUL15
[1] "16"

$CMPORTNE09UW
[1] "8"

$CMROUGH09UW
 [1] "1"  "5"  "13" "25" "30" "33" "35" "44" "46" "51" "64" "70" "74" "78" "86" "87" "89" "93"

$CMRUB96UW
[1] "31" "71" "83"

$CMRUDY10
[1] "61" "92"

$CMRUS09UW
[1] "48"

$CMRUS93UW
[1] "74" "79"

$CMRUSSI07UW
[1] "48"

$CMRUSSI09UW
[1] "50" "63" "69"

$CMSANC96UW
[1] "6"  "89"

$CMSITKI09UW
[1] "87" "95"

$CMSPIRU09UW
[1] "30" "38" "42" "44" "58" "93"

$CMSTE09UW
[1] "69"

$CMSTEPR09UW
[1] "95"

$CMSTU09UW
 [1] "64"  "105" "106" "107" "108" "109" "110" "111" "112" "113" "114" "115" "116" "117" "118" "119" "120" "121" "122" "123"
[21] "124" "125"

$CMSTU92
 [1] "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10" "11" "12" "13" "14" "15" "16" "17" "18" "19" "20" "21" "22" "23" "24"
[25] "25" "26" "27" "28" "29" "30" "31" "32" "33" "34" "35" "36" "37" "38" "39" "40" "41" "42" "43" "44" "45" "46" "47" "48"
[49] "49" "50" "51" "52" "53" "54" "55" "56" "57" "58" "59" "60" "61" "62" "63" "64" "65" "66" "67" "68" "69" "70" "71"

$CMUGAN09UW
 [1] "11" "16" "30" "43" "45" "49" "51" "63" "76" "86" "87" "88" "91" "95"

$CMVOL96UW
[1] "27"

$CMWESA93UW
[1] "6"

$CMWESE93UW
[1] "6"

$CMWESF93UW
[1] "44" "45" "46" "47" "48" "68" "86"

$CMWESG93UW
[1] "4"

$CMWESH93UW
[1] "34" "73"

$CMWESI93UW
[1] "55"

$CMWESK93UW
[1] "30"

$CMWESL93UW
[1] "32"

$CMWESN93UW
[1] "23" "36"

$CMWESP93UW
[1] "17" "37" "46"

$CMWKIL09UW
[1] "34" "42" "59" "70"

$CMZAC92UW
[1] "68" "71" "74"

$CMZACR09UW
[1] "37"

sample_size_qa <- sample_size_qa %>% 
  mutate(missing = genotyped - sapply(silly, function(x) get(paste0(x, ".gcl"))$n))

Duplicate

Remove duplicate individuals within the same collection. Typically we specify duplicates as a pair of individuals that share >=95% of genotypes. Once a pair of duplicates is identified, we keep the individual with the most genotypes and remove the other.

duplicate_check_95 <- CheckDupWithinSilly.GCL(sillyvec = collections_91, loci = loci96, quantile = NULL, minproportion = 0.95)
duplicate_summary <- sapply(collections_91, function(x) duplicate_check_95[[x]]$report, simplify = FALSE)
duplicate_remove <- RemoveDups.GCL(duplicate_check_95)
save_objects(c("duplicate_summary", "duplicate_remove"), "../2018/Objects/")
# show individuals removed
duplicate_remove[duplicate_remove != "Nothing Removed"]

$`CMBARLR15E`
[1] "87"  "21"  "100" "13" 

$CMBEARBC09UW
[1] "38"

$CMBSUK09UW
[1] "36"

$CMCHI09UW
[1] "19" "34" "44"

$CMCHI96UW
[1] "7"

$CMCHIGK09UW
[1] "48"

$CMCOXC09UW
[1] "23"

$CMEAGLH09UW
[1] "30"

$CMIVAN09UW
[1] "5"

$CMKARL09UW
[1] "80"

$CMKIAL09UW
[1] "11"

$CMKIAV17
[1] "30"

$CMKIZH09UW
[1] "58"

$CMOLDS10
[1] "38" "80"

$CMRUSSI07UW
[1] "14"

$CMSANC96UW
[1] "49" "36"

$CMSTE09UW
[1] "32"

$CMWESN93UW
[1] "10"

$CMZAC92UW
[1] "13"

sample_size_qa <- sample_size_qa %>% 
  mutate(duplicate = genotyped - missing - sapply(silly, function(x) get(paste0(x, ".gcl"))$n))

Final

How many fish did we end up with in our final baseline? Save final genotypes.

(sample_size_qa <- sample_size_qa %>% 
  mutate(final = sapply(silly, function(x) get(paste0(x, ".gcl"))$n)))

write_csv(sample_size_qa, "../2018/Tables/sample_size_qa.csv")
dir.create("../2018/Genotypes/original_split_postQA/")
save_sillys(collections_91, "../2018/Genotypes/original_split_postQA/")

What is smallest collection?

Whoa, we should totally junk some of these collections given their low sample size (or complete lack of fish!). We’ll keep anything with at least 40 fish.

(collections_88 <- sample_size_qa %>% 
   filter(final > 40) %>% 
   pull("silly"))

 [1] "CHAME92UW"    "CHBAL92UW"    "CHDOG92UW"    "CHKIZ92UW"    "CHLIJ92UW"    "CHUGA92UW"    "CMAMBM09UW"   "CMBARL09UW"  
 [9] "CMBARL17"     "CMBARLR15E"   "CMBARLR15L"   "CMBARLR15M"   "CMBEARBC09UW" "CMBEL92UW"    "CMBIGRI09UW"  "CMBSU92UW"   
[17] "CMBSUK09UW"   "CMCAN09UW"    "CMCAN92UW"    "CMCHI09UW"    "CMCHI96UW"    "CMCHIGK09UW"  "CMCOL96UW"    "CMCOXC09UW"  
[25] "CMDEAD09UW"   "CMDEL96UW"    "CMDOGSC10"    "CMDRYBR09UW"  "CMEAGLH09UW"  "CMFOST09UW"   "CMGULLC09UW"  "CMHIDDB15"   
[33] "CMIVAN09UW"   "CMKAIUG15"    "CMKARL09UW"   "CMKIAL09UW"   "CMKIAV17"     "CMKITB09UW"   "CMKIZH09UW"   "CMKUJUNF09UW"
[41] "CMLIJ09UW"    "CMMIDWAY15"   "CMNAKIL08UW"  "CMNATA17"     "CMNECR08UW"   "CMOCEB09UW"   "CMOLDS10"     "CMPASS09UW"  
[49] "CMPAUL15"     "CMPORTC08UW"  "CMPORTNE09UW" "CMROUGH09UW"  "CMRUB96UW"    "CMRUDY10"     "CMRUS09UW"    "CMRUS92UW"   
[57] "CMRUS93UW"    "CMRUSSI07UW"  "CMRUSSI09UW"  "CMSANC09UW"   "CMSANC96UW"   "CMSITKI09UW"  "CMSMOKH10"    "CMSPIRU09UW" 
[65] "CMSTE09UW"    "CMSTE92UW"    "CMSTEPR09UW"  "CMSTU09UW"    "CMUGAN09UW"   "CMVOL09UW"    "CMVOL92UW"    "CMVOL96UW"   
[73] "CMWESA93UW"   "CMWESD93UW"   "CMWESE93UW"   "CMWESF93UW"   "CMWESG93UW"   "CMWESH93UW"   "CMWESI93UW"   "CMWESJ93UW"  
[81] "CMWESK93UW"   "CMWESL93UW"   "CMWESM93UW"   "CMWESN93UW"   "CMWESP93UW"   "CMWKIL09UW"   "CMZAC92UW"    "CMZACR09UW"  

save_objects("collections_88", "../2018/Objects/")

Check HWE - collections

Write out a Genepop file to check HWE within collections.

dir.create("../2018/Genepop")
gcl2Genepop.GCL(sillyvec = collections_88, loci = loci96, path = "../2018/Genepop/collections_88.gen", VialNums = TRUE)

NULL

Calculate HWE using the genepop package link.

test_HW(inputFile = "../2018/Genepop/collections_88.gen", outputFile = "../2018/Genepop/collections_88.P.txt")

Dangerous input file name extension:
   the genepop code expects either extension '.txt' or no extension.

Results are stored in file ../2018/Genepop/collections_88.P.txt 

Summary table

Read in results and save summary HWE p-values.

hwe_collections_88 <- ReadGenepopHWE.GCL(file = "../2018/Genepop/collections_88.txt.P", sillyvec = collections_88)

Read 18222 items

All p-values are pulled directly from the Genepop *.P file if calculated via Exact test.
If calculated via MCMC then over all loci/pops are derived in this function, correcting p = 0 for the number of batches * iterations!
NA in 'smmry' dataframe means that locus in either 1) monomorphic (or very low MAF) or 2) not diploid (haploid mtSNP).
Haploid loci are not included in the 'HWdata' data.frame.NULL

(hwe_summary_collections_88 <- hwe_collections_88$SummaryPValues %>% 
  as_tibble(rownames = "locus"))

write_csv(x = hwe_summary_collections_88, "../2018/Tables/hwe_summary_collections_88.csv")

Collections out of HWP

What do overall-loci p-values look like for collections?

All look pretty good with the potential exception of CMBARLR15L.

Loci out of HWP

What do overall-loci p-values look like for loci?

There are a couple of loci that are suspect.

HWP plots

Per Waples 2014, we will plot a histogram of HWP p-values to look for collections/loci that do not conform to HWP. We will limit our search to the specific collection and loci that were suspicious.

Collection

Plot p-values per locus for CMBARLR15L

hwe_summary_collections_88 %>% 
  gather(silly, p, -locus) %>%
  filter(silly == "CMBARLR15L") %>% 
  filter(locus != "Overall Loci") %>% 
  ggplot(aes(p)) + 
  geom_histogram(binwidth = 0.05) +
  geom_hline(yintercept = 5, colour = "red") +
  ggtitle("HWP p-values by locus for CMBARLR15L")

Looks fine, no major cause for concern, now on to plotting for loci.

Loci

Plot p-values per collection for loci with overall pops p-value < 0.2.

# filter for loci with overall pops < 0.2
hwe_loci_plot <- hwe_summary_collections_88 %>% 
  gather(silly, p, -locus) %>%
  filter(silly == "Overall Pops" & p < 0.2) %>% 
  pull(locus)
# plot
hwe_summary_collections_88 %>% 
  gather(silly, p, -locus) %>%
  filter(locus %in% hwe_loci_plot) %>% 
  filter(silly != "Overall Pops") %>%
  ggplot(aes(p)) + 
  geom_histogram(binwidth = 0.05) +
  geom_hline(yintercept = length(collections_88) / 20, colour = "red") +
  ggtitle("HWP p-values by collection") +
  facet_grid( ~ locus)

Oke_U2041-84 and Oke_U506-110 are a bit suspect, but we’ll let them go since their overall p-values were > 0.05.

Takeaway

All collections and loci were retained.

Pooling temporal collections

Now I’ll follow Andy’s pooling of temporal collections in to populations.

Non-Barling

Andy’s

Create a list for pooling (excluding Barling). Note some are different than Andy’s (i.e. Stepovak)

temporal_pooling <- list("Amber Bay" = c("CMWESN93UW", "CMAMBM09UW"),
                         "Balboa Bay" = c("CHBAL92UW", "CMFOST09UW"),
                         "Hallo Bay" = c("CMWESH93UW", "CMBIGRI09UW"),
                         "Bear Bay" = c("CMWESD93UW", "CMBEARBC09UW"),
                         "Big Sukhoi" = c("CMBSU92UW", "CMBSUK09UW"),
                         "Canoe Bay" = c("CMCAN92UW", "CMCAN09UW"),
                         "Chichagof Bay" = c("CMCHI96UW", "CMCHI09UW"),
                         "Chigniagak Bay" = c("CMWESI93UW", "CMCHIGK09UW"),
                         "Gull Cape Creek" = c("CMWESP93UW", "CMGULLC09UW"),
                         "Ivanof River" = c("CMWESL93UW", "CMIVAN09UW"),
                         "Kialagvik Creek" = c("CMWESF93UW", "CMKIAL09UW"),
                         "Kitoi Bay" = c("CMWESA93UW", "CMKITB09UW"),
                         "Kizhuyak River" = c("CHKIZ92UW", "CMKIZH09UW"),
                         "Kujulik Bay" = c("CMWESK93UW", "CMKUJUNF09UW"),
                         "Little John Lagoon" = c("CHLIJ92UW", "CMLIJ09UW"),
                         "Portage Creek" = c("CMWESJ93UW", "CMPORTC08UW"),
                         "Russell Creek" = c("CMRUS92UW", "CMRUS93UW", "CMRUS09UW"),
                         "Russian River" = c("CMRUSSI07UW", "CMRUSSI09UW"),
                         "Sandy Cove" = c("CMSANC96UW", "CMSANC09UW"),
                         "Stepovak Bay" = c("CMSTE92UW", "CMWESM93UW", "CMSTE09UW"),
                         "Uganik River" = c("CHUGA92UW", "CMUGAN09UW"),
                         "Volcano Bay" = c("CMVOL92UW", "CMVOL96UW", "CMVOL09UW")
)
save_objects("temporal_pooling", "../2018/Objects/")
# verify that all sillys are in collections_88
setdiff(unlist(temporal_pooling), collections_88)

character(0)

Now that we have our list of putative populations to pool, let’s see how they do in Fisher’s test.

collections_88_allele_freq <- FreqPop.GCL(sillyvec = collections_88, loci = loci96)
temporal_pooling_results <- FishersTest.GCL(freq = collections_88_allele_freq, loci = loci96, tests = temporal_pooling)
save_objects(c("collections_88_allele_freq", "temporal_pooling_results"), "../2018/Objects/")
# overall loci p-values
temporal_pooling_results$OverallResults %>% 
  as_tibble(rownames = "pop") %>% 
  arrange(overall)

WASSIP

I went and double checked the WASSIP chum baseline SP 12-26, Table 2. Should have done this first.

• The 1992 and 1996 Volcano Bay collections got dropped (-CMVOL92UW, -CMVOL96UW)
• Pool Stepovak Bay 1992 and 2009, keep Stepovak River 2009 separate, drop 1993 Stepovak River (-CMWESM93UW)
• Drop Amber Bay 1993 (-CMWESN93UW)
• Drop Hallo Bay 1993 (-CMWESH93UW)

temporal_pooling_WASSIP <- list("Balboa Bay" = c("CHBAL92UW", "CMFOST09UW"),
                                "Bear Bay" = c("CMWESD93UW", "CMBEARBC09UW"),
                                "Big Sukhoi" = c("CMBSU92UW", "CMBSUK09UW"),
                                "Canoe Bay" = c("CMCAN92UW", "CMCAN09UW"),
                                "Chichagof Bay" = c("CMCHI96UW", "CMCHI09UW"),
                                "Chigniagak Bay" = c("CMWESI93UW", "CMCHIGK09UW"),
                                "Gull Cape Creek" = c("CMWESP93UW", "CMGULLC09UW"),
                                "Ivanof River" = c("CMWESL93UW", "CMIVAN09UW"),
                                "Kialagvik Creek" = c("CMWESF93UW", "CMKIAL09UW"),
                                "Kitoi Bay" = c("CMWESA93UW", "CMKITB09UW"),
                                "Kizhuyak River" = c("CHKIZ92UW", "CMKIZH09UW"),
                                "Kujulik Bay" = c("CMWESK93UW", "CMKUJUNF09UW"),
                                "Little John Lagoon" = c("CHLIJ92UW", "CMLIJ09UW"),
                                "Portage Creek" = c("CMWESJ93UW", "CMPORTC08UW"),
                                "Russell Creek" = c("CMRUS92UW", "CMRUS93UW", "CMRUS09UW"),
                                "Russian River" = c("CMRUSSI07UW", "CMRUSSI09UW"),
                                "Sandy Cove" = c("CMSANC96UW", "CMSANC09UW"),
                                "Stepovak Bay" = c("CMSTE92UW", "CMSTE09UW"),
                                "Uganik River" = c("CHUGA92UW", "CMUGAN09UW")
)
save_objects("temporal_pooling_WASSIP", "../2018/Objects/")
# verify that all sillys are in collections_88
setdiff(unlist(temporal_pooling_WASSIP), collections_88)

character(0)

Now that we have our list of WASSIP populations to pool, let’s see how they do in Fisher’s test.

temporal_pooling_WASSIP_results <- FishersTest.GCL(freq = collections_88_allele_freq, loci = loci96, tests = temporal_pooling_WASSIP)
save_objects("temporal_pooling_WASSIP_results", "../2018/Objects/")
# overall loci p-values
temporal_pooling_WASSIP_results$OverallResults %>% 
  as_tibble(rownames = "pop") %>% 
  arrange(overall)

Looks good enough to me, also coincides with WASSIP, so that is good news. Now on to check Barling.

Barling

Test all possible pairwise combinations of Barling samples

barling_pooling <- combn(x = grep(pattern = "BARL", x = collections_88, value = TRUE), m = 2, simplify = FALSE)
save_objects("temporal_pooling_WASSIP", "../2018/Objects/")
barling_pooling_results <- FishersTest.GCL(freq = collections_88_allele_freq, loci = loci96, tests = barling_pooling)
save_objects("barling_pooling_results", "../2018/Objects/")
# overall loci p-values
barling_pooling_results$OverallResults %>% 
  as_tibble(rownames = "pop") %>% 
  arrange(overall)

While it certainly does seem like we “could” pool some of the Barling collections per our WASSIP rules, I think it is best to leave them all separate for now for plotting purposes.

Final pooling

Accept all WASSIP pooling and move on.

invisible(
  lapply(temporal_pooling_WASSIP, function(pop) {
    PoolCollections.GCL(collections = pop, loci = loci96, IDs = NULL, newname = paste(pop, collapse = "."))
  } )
)

Ordered populations

Need new sillyvec with pooled populations.

collections_post_pooling <- c(
  setdiff(setdiff(collections_88, 
          unlist(temporal_pooling_WASSIP)),  # get single collections from non-WASSIP pooling
  c("CMVOL92UW", "CMVOL96UW", "CMWESM93UW", "CMWESN93UW", "CMWESH93UW")), # drop these old collections that didn't pool
  sapply(temporal_pooling_WASSIP, function(pop) {
    paste(pop, collapse = ".")  # add pooled WASSIP
  })
)
dir.create("../2018/Genotypes/original_split_postQA_postpooling/")

'..\2018\Genotypes\original_split_postQA_postpooling' already exists

save_sillys(collections_post_pooling, "../2018/Genotypes/original_split_postQA_postpooling/")

Error in save_sillys(collections_post_pooling, "../2018/Genotypes/original_split_postQA_postpooling/") : 
  These sillys:
CHAME92UW
CHDOG92UW
CMAMBM09UW
CMBARL09UW
CMBARL17
CMBARLR15E
CMBARLR15L
CMBARLR15M
CMBEL92UW
CMBIGRI09UW
CMCOL96UW
CMCOXC09UW
CMDEAD09UW
CMDEL96UW
CMDOGSC10
CMDRYBR09UW
CMEAGLH09UW
CMHIDDB15
CMKAIUG15
CMKARL09UW
CMKIAV17
CMMIDWAY15
CMNAKIL08UW
CMNATA17
CMNECR08UW
CMOCEB09UW
CMOLDS10
CMPASS09UW
CMPAUL15
CMPORTNE09UW
CMROUGH09UW
CMRUB96UW
CMRUDY10
CMSITKI09UW
CMSMOKH10
CMSPIRU09UW
CMSTEPR09UW
CMSTU09UW
CMVOL09UW
CMWESE93UW
CMWESG93UW
CMWKIL09UW
CMZAC92UW
CMZACR09UW
CHBAL92UW.CMFOST09UW
CMWESD93UW.CMBEARBC09UW
CMBSU92UW.CMBSUK09UW
CMCAN92UW.CMCAN09UW
CMCHI96UW.CMCHI09UW
CMWESI93UW.CMCHIGK09UW
CMWESP93UW.CMGULLC09UW
CMWESL93UW.CMIVAN09UW
CMWESF93UW.CMKIAL09UW
CMWESA93UW.CMKITB09UW
CHKIZ92UW.CMKIZH09UW
CMWESK93UW.CMKUJUNF09UW
CHLIJ92UW.CMLIJ09UW
CMWESJ93UW.CMPORTC08UW
CMRUS92UW.CMRUS93UW.CMRUS09UW
CMRUSSI07UW.CMRUSSI09UW
CMSANC96UW.CMSANC09UW
CMSTE92UW.CMSTE09UW
CHUGA92UW.CMUGAN09UW
are not in your workspace, hoser!!!

Now need to join up these post pooling collections (populations) with the paired OceanAK data from collections_89_info and get in geographical order.

Add new grouping variable for:
* South Peninsula * Chignik * Kodiak Mainland * Kodiak * Sturgeon/Kitoi * Barling River * 2017 Collections

populations_63_info <- tibble(pop = collections_post_pooling) %>% 
  mutate(n = sapply(collections_post_pooling, function(x) {get(paste0(x, ".gcl"))$n} )) %>%  # get final sample sizes
  mutate(pop_no_UW = str_remove_all(pop, "UW")) %>%  # remove UW extension
  mutate(pop_no_UW = dplyr::recode(pop_no_UW, !!!list("CMBARLR15E" = "CMBARLR15", "CMBARLR15M" = "CMBARLR15", "CMBARLR15L" = "CMBARLR15"))) %>%  # remove E, M, L extension for Barling 2015 for joining
  separate(pop_no_UW, c("col1", "col2", "col3"), sep = "\\.", remove = TRUE) %>%  # separate sillys to join dates
  left_join(select(collections_89_info, `Silly Code`, `Collection Date`), by = c("col1" = "Silly Code")) %>%  # get date 1
  left_join(select(collections_89_info, `Silly Code`, `Collection Date`), by = c("col2" = "Silly Code")) %>%  # get date 2
  left_join(select(collections_89_info, `Silly Code`, `Collection Date`), by = c("col3" = "Silly Code")) %>%  # get date 3
  rename("Collection Date 1" = "Collection Date.x", "Collection Date 2" = "Collection Date.y", "Collection Date 3" = "Collection Date") %>%  # rename dates
  mutate(`Collection Date 1` = str_replace(`Collection Date 1`, "1/1/1992", "1992")) %>%  # CMZAC09 has bogus date
  mutate(col = ifelse(is.na(col3), ifelse(is.na(col2), col1, col2), col3)) %>%  # get latest collection from each pop
  left_join(select(collections_89_info, `Silly Code`, Quadrant, color, Location, Latitude, Longitude), by = c("col" = "Silly Code")) %>%  # join OceanAK based on latest collection from each pop
  mutate(Group = case_when(Quadrant == "Peninsula - South" ~ "South Peninsula",
                           Quadrant == "Chignik" ~ "Chignik",
                           Quadrant == "Kodiak Mainland" ~ "Kodiak Mainland",
                           pop %in% c("CMPAUL15", "CMSTU09UW", "CMWESA93UW.CMKITB09UW") ~ "Sturgeon/Kitoi Bay",
                           pop %in% c("CMBARL09UW", "CMBARLR15E", "CMBARLR15M", "CMBARLR15L") ~ "Barling Bay",
                           pop %in% c("CMNATA17", "CMKIAV17", "CMBARL17") ~ "2017 Collections",
                           TRUE ~ "Kodiak")) %>%  # create group for groupvec
  dplyr::select(pop, Location, `Collection Date 1`, `Collection Date 2`, `Collection Date 3`, Group, Quadrant, n, Latitude, Longitude, color) %>%  # keep columns
  separate(Location, c("Location", "trash"), sep = " - ") %>%  # standardize Locations part 1
  separate(Location, c("Location", "more_trash"), sep = " \\(") %>%  # standardize Locations part 2
  dplyr::select(-trash, -more_trash) %>%  # drop bogus location info
  mutate(`Collection Date 1` = case_when(pop == "CMBARLR15M" ~ "8/3/2015",
                                         pop == "CMBARLR15L" ~ "9/11/2015",
                                         TRUE ~ `Collection Date 1`)) %>%  # fix Barling 2015 dates
  mutate(Month = month(mdy(`Collection Date 1`), label = TRUE)) %>% #  extract month
  dplyr::select(pop, Location, `Collection Date 1`, `Collection Date 2`, `Collection Date 3`, Month, Group, Quadrant, n, Latitude, Longitude, color)  # keep columns

Create a leaflet map to determine map order for populations.

populations_63_info <- populations_63_info %>% 
  mutate(color = case_when(Group == "South Peninsula" ~ "red",
                           Group == "Chignik" ~ "darkred",
                           Group == "Kodiak Mainland" ~ "orange",
                           Group == "Kodiak" ~ "green",
                           Group == "Sturgeon/Kitoi Bay" ~ "blue",
                           Group == "Barling Bay" ~ "purple",
                           Group == "2017 Collections" ~ "darkpurple"))
colors_7 <- c("red", "darkred", "orange", "blue", "green", "purple", "violet")
save_objects("colors_7", "../2018/Objects/")
icons <- awesomeIcons(icon = 'circle', iconColor = 'black', library = 'ion', markerColor = populations_63_info$color)
populations_63_info %>% 
  leaflet(width = "100%") %>% 
  addTiles() %>% 
  addAwesomeMarkers(~ Longitude, ~ Latitude, icon = icons, label = ~ pop, labelOptions = labelOptions(noHide = TRUE, direction = "auto"))

Save geography ordered population vector.

populations_63 <- read_csv("../2018/populations_63.csv")

Parsed with column specification:
cols(
  pop_no = col_integer(),
  map_no = col_integer(),
  pop = col_character()
)

populations_63_info <- populations_63_info %>% 
  left_join(populations_63, by = "pop") %>% 
  select(pop_no, map_no, pop, Location, `Collection Date 1`, `Collection Date 2`, `Collection Date 3`, Month, Group, Quadrant, n, Latitude, Longitude, color) %>% 
  arrange(pop_no)
write_csv(populations_63_info, "../2018/Tables/populations_63_info.csv")
populations_63 <- populations_63$pop
groups_7 <- unique(populations_63_info$Group)
groupvec_7 <- match(populations_63_info$Group, groups_7)
save_objects(c("populations_63", "populations_63_info", "groups_7", "groupvec_7"), "../2018/Objects/")

Final map

Make map with the correct map numbers to double check.

icons <- awesomeIcons(icon = 'egg', iconColor = 'transparent', library = 'ion', markerColor = populations_63_info$color)
populations_63_info %>% 
  leaflet(width = "100%") %>% 
  addTiles() %>% 
  addAwesomeMarkers(~ Longitude, ~ Latitude, icon = icons, label = ~ as.character(map_no), popup = ~ Location, labelOptions = labelOptions(noHide = TRUE, textOnly = TRUE, direction = "top"))

Check HWE - populations

Write out a Genepop file to check HWE within populations.

gcl2Genepop.GCL(sillyvec = populations_63, loci = loci96, path = "../2018/Genepop/populations_63.gen", VialNums = TRUE)

NULL

Calculate HWE using the genepop package link.

test_HW(inputFile = "../2018/Genepop/populations_63.gen", outputFile = "../2018/Genepop/populations_63.txt.P")

Dangerous input file name extension:
   the genepop code expects either extension '.txt' or no extension.

Results are stored in file ../2018/Genepop/populations_63.txt.P 

Summary table

Read in results and save summary HWE p-values.

hwe_populations_63 <- ReadGenepopHWE.GCL(file = "../2018/Genepop/populations_63.txt.P", sillyvec = populations_63)

Loading required package: reshape2

Attaching package: 㤼㸱reshape2㤼㸲

The following object is masked from 㤼㸱package:tidyr㤼㸲:

    smiths

Read 13322 items

All p-values are pulled directly from the Genepop *.P file if calculated via Exact test.
If calculated via MCMC then over all loci/pops are derived in this function, correcting p = 0 for the number of batches * iterations!
NA in 'smmry' dataframe means that locus in either 1) monomorphic (or very low MAF) or 2) not diploid (haploid mtSNP).
Haploid loci are not included in the 'HWdata' data.frame.NULL

(hwe_summary_populations_63 <- hwe_populations_63$SummaryPValues %>% 
  as_tibble(rownames = "locus"))

write_csv(x = hwe_summary_populations_63, "../2018/Tables/hwe_summary_populations_63.csv")

Populations out of HWP

What do overall-loci p-values look like for collections?

hwe_summary_populations_63 %>% 
  gather(silly, p, -locus) %>%
  filter(locus == "Overall Loci") %>% 
  arrange(p)

Looks pretty good, nothing major to worry aobut. On to Fst!

Combine loci

Planning to use WASSIP locus information (ADF&G SP12-26):

• Combing 3 haploid loci (mitochondrial)
  • Oke_Cr30
  • Oke_Cr386
  • Oke_ND3-69
• Combine 2 linked diploid loci (nuclear)
  • Oke_U1021-102
  • Oke_U1022-139
• Drop 2 diploid loci (nuclear) due to linkage
  • Oke_gdh1-62
  • Oke_pgap-92

mito_combine <- names(which(LocusControl$ploidy == 1))
nuclear_combine <- c("Oke_U1021-102", "Oke_U1022-139")
nuclear_drop <- c("Oke_gdh1-62", "Oke_pgap-92")
CombineLoci.GCL(sillyvec = populations_63, markerset = nuclear_combine, update = TRUE, delim = ".")

NULL

CombineLoci.GCL(sillyvec = populations_63, markerset = mito_combine, update = TRUE, delim = ".")

NULL

dir.create("../2018/Genotypes/original_split_postQA_postpooling_combineloci/")
save_sillys(populations_63, "../2018/Genotypes/original_split_postQA_postpooling_combineloci/")
loci91 <- c(
  setdiff(loci96, 
          c(nuclear_drop, mito_combine, nuclear_combine)),  # drop linked markers and single combined
  paste(nuclear_combine, collapse = "."), paste(mito_combine, collapse = "."))  # add combined

Allele frequency plots

Calculate allele frequencies.

dir.create("../2018/FreqPlots")
allele_freq <- FreqFisPlot4SNPs.GCL(sillyvec = populations_63, loci = loci96, groupvec = groupvec_7, groupcol = colors_7, file = "../2018/FreqPlots/allele_freq_populations_63_loci96.pdf")

Fst

Calculate

Now that we have our final post-QA, post-pooling genotypes, time to calculate Fst between all of our populations.

dir.create("../2018/FSTAT")
Fst_tree_populations_63_loci91 <- PairwiseFstTree.GCL(sillyvec = populations_63, loci = loci91, dir = "../2018/FSTAT", nboots = 1000, ncores = 4, returnbootstrapFst = FALSE)

Heatmap of Fst

pairwise_fst <- Fst_tree_populations_63_loci91$PairwiseFst
write_csv(as_data_frame(pairwise_fst, rownames = "pop"), "../2018/Tables/pairwise_fst_populations_63.csv")
save_objects("pairwise_fst", "../2018/Objects/")
get_lower_tri <- function(mat) {
  mat[upper.tri(mat)] <- NA
  return(mat)
}
get_lower_tri(pairwise_fst) %>% 
  as_tibble(rownames = "pop1") %>%
  mutate(pop1 = factor(pop1, populations_63)) %>% 
  gather(pop2, fst, -pop1, na.rm = TRUE, factor_key = TRUE) %>% 
  ggplot(aes(x = pop1, y = pop2, fill = fst)) +
  geom_tile() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  ggtitle("Pairwise Fst - 63 Populations, 91 Loci")

MDS

All populations

Make an MDS of the Fst data for comparison to Andy’s work back in 2015.

dist = pairwise_fst
popvec = 1:63
groupsvec = 1:7
colvec = match(colors_7[groupvec_7], colors())
groups = groups_7[groupsvec]
cols = colors_7[groupsvec]
main = ""
labels = TRUE
locNames = 1:63
axes = TRUE
box = TRUE
adj = c(1.5, 1.5)
cex = 1
size = 1.3
#~~~~~~~~
names <- if(labels==F){NULL} else{if(locNames==F){popvec} else{locNames[popvec]}}
xx=dist[popvec,popvec]
x=as.vector(cmdscale(xx,k=3)[,1])
y=as.vector(cmdscale(xx,k=3)[,2])
z=-as.vector(cmdscale(xx,k=3)[,3])
plot3d(x,y,z+abs(range(z)[1]),aspect=F,col=colors()[colvec[popvec]],size=size,type='s',main=main,box=box,axes=axes,top=T,cex=1,xlab='',ylab='',zlab='',xlen=0,ylen=0,zlen=0)
plot3d(x,y,z+abs(range(z)[1]),aspect=F,col="black",size=3,type='h',box=F,axes=F,top=T,add=T,xlab='',ylab='',zlab='',xlen=0,ylen=0,zlen=0)
if(labels==T){texts3d(x,y,z+abs(range(z)[1]),adj=adj,text=names,font=1,cex=cex,add=T,top=T,axes=F,xlab='',ylab='',zlab='')}
par3d(windowRect = c(0, 0, 2000, 2000))
legend3d("topright", legend = groups, pch = 16, col = cols, inset = c(0), cex = 5)
rglwidget()

Kodiak populations

Make an MDS of the Fst data for just Kodiak populations

dist = pairwise_fst
popvec = 32:63
colvec = match(colors_7[groupvec_7], colors())
groups = groups_7[4:7]
cols = colors_7[4:7]
main = ""
labels = TRUE
locNames = 1:63
axes = TRUE
box = TRUE
adj = c(1.5, 1.5)
cex = 1
size = 1.3
#~~~~~~~~
names <- if(labels==F){NULL} else{if(locNames==F){popvec} else{locNames[popvec]}}
xx=dist[popvec,popvec]
x=as.vector(cmdscale(xx,k=3)[,1])
y=as.vector(cmdscale(xx,k=3)[,2])
z=-as.vector(cmdscale(xx,k=3)[,3])
plot3d(x,y,z+abs(range(z)[1]),aspect=F,col=colors()[colvec[popvec]],size=size,type='s',main=main,box=box,axes=axes,top=T,cex=1,xlab='',ylab='',zlab='',xlen=0,ylen=0,zlen=0)
plot3d(x,y,z+abs(range(z)[1]),aspect=F,col="black",size=3,type='h',box=F,axes=F,top=T,add=T,xlab='',ylab='',zlab='',xlen=0,ylen=0,zlen=0)
if(labels==T){texts3d(x,y,z+abs(range(z)[1]),adj=adj,text=names,font=1,cex=cex,add=T,top=T,axes=F,xlab='',ylab='',zlab='')}
par3d(windowRect = c(0, 0, 2000, 2000))
legend3d("topright", legend = groups, pch = 16, col = cols, inset = c(0), cex = 5)
rglwidget()

Kodiak populations, no outliers

Remove Sturgeon/Kito Bay populations, Karluk Lagoon, Gull Cape Lagoon, Big Sukhoi, and Sitkinak Island as outliers.

dist = pairwise_fst
popvec = populations_63_info %>% 
  mutate(pop_no = 1:63) %>% 
  filter(Quadrant == "Kodiak/Afognak Islands") %>% 
  filter(Group != "Sturgeon/Kitoi Bay") %>% 
  filter(!Location %in% c("Karluk Lagoon", "Gull Cape Lagoon", "Big Sukhoi", "Sitkinak Island")) %>% 
  pull(pop_no)
groupsvec = 5:7
colvec = match(colors_7[groupvec_7], colors())
groups = groups_7[groupsvec]
cols = colors_7[groupsvec]
main = ""
labels = TRUE
locNames = 1:63
axes = TRUE
box = TRUE
adj = c(1.5, 1.5)
cex = 1
size = 1.3
#~~~~~~~~
names <- if(labels==F){NULL} else{if(locNames==F){popvec} else{locNames[popvec]}}
xx=dist[popvec,popvec]
x=as.vector(cmdscale(xx,k=3)[,1])
y=as.vector(cmdscale(xx,k=3)[,2])
z=-as.vector(cmdscale(xx,k=3)[,3])
plot3d(x,y,z+abs(range(z)[1]),aspect=F,col=colors()[colvec[popvec]],size=size,type='s',main=main,box=box,axes=axes,top=T,cex=1,xlab='',ylab='',zlab='',xlen=0,ylen=0,zlen=0)
plot3d(x,y,z+abs(range(z)[1]),aspect=F,col="black",size=3,type='h',box=F,axes=F,top=T,add=T,xlab='',ylab='',zlab='',xlen=0,ylen=0,zlen=0)
if(labels==T){texts3d(x,y,z+abs(range(z)[1]),adj=adj,text=names,font=1,cex=cex,add=T,top=T,axes=F,xlab='',ylab='',zlab='')}
par3d(windowRect = c(0, 0, 2000, 2000))
legend3d("topright", legend = groups, pch = 16, col = cols, inset = c(0), cex = 5)
rglwidget()

Summary

Added 3 new collections:

1. CMBARL17 - Barling Bay Creek - 9/1/2017
2. CMKIAV17 - Kiavak Portage - 9/1/2017
3. CMNATA17 - Natalia Bay Creek - 9/1/2017

Of note, the 2017 collection from Barling Bay Creek came from higher up in the drainaage (Birch Foster, personal communication), but still looks very similar to the 2015 September collection. Kiavak Portage is closely related to Barling, and Natalia was bit further away from the main cluster.

This update does not change any departmental recommendations regarding the use of local broodsource (Barling Bay Creek) for a chum hatchery in Old Harbor.